Soft modes and elasticity of nearly isostatic lattices: randomness and dissipation 
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The square lattice with nearest neighbor central-force springs is isostatic and does not support 
shear. Using the Coherent Potential Approximation (CPA), we study how the random addition, with 
probability V = (z— 4) /4 (z = average number of contacts), of next-nearest-neighbor (NNN) springs 
restores rigidity and affects phonon structure. The CPA effective NNN spring constant « m (w), 
equivalent to the complex shear modulus G(u), obeys the scaling relation, k m (ui) — K m h(to /ui*), at 
small V, where K m = k' m (0) ~ V 2 and uj* ~ V, implying nonafhne elastic response at small V and 
the breakdown of plane- wave states beyond the loffe-Regel limit at uj ~ uj* . We identify a divergent 
length I* ~ V~ , and we relate these results to jamming. 

PACS numbers: 61.43.-j, 62.20.de, 46.65.+g, 05.70.Jk 



Isostatic lattices [3-0| are systems at the onset of me- 
chanical stability in which the average number of con- 
tacts z per particle in d-dimensions is equal to z c = 2d. 
A lattice with TV particles and N c two-particle contacts 
has No = dN — N c zero modes. An infinite isostatic 
lattice is one in which N c = Nz c /2, and the fraction of 
zero modes vanishes. Because particles at the boundary 
have fewer contacts than those in the bulk, the num- 
ber of zero modes in a finite isostatic lattice is subex- 
tensive (N ~ N( d - x )/ d ) and proportional to the area of 
the system boundary. As a result, the phonon spectrum 
of isostatic lattices is one-dimensional in nature. These 
properties underly the elastic and vibrational properties 
of a variety of systems including network glasses 0, [B| , 
rigidity percolation 0,0 j /3-cristobalite granular me- 
dia [H, , and networks of semi-flexible polymers [n} . 
Isostatic lattices include d-dimensional hypercubic lat- 
tices and the 2d kagome, the 3d pyrochlore lattice, and 
their d-dimensional generalizations |12[ , all with central- 
force springs with spring constant k connecting nearest 
neighbor (NN) sites. They also include randomly packed 
spheres at the jamming transition [uj - flo| . 

As in critical phenomena at "standard" phase transi- 
tions, the approach to the critical isostatic state, which 
this paper explores, is characterized by diverging length 
and time scales and by scaling behavior. Lattices can 
be moved off isostaticity in various ways, including (1) 
introducing springs with a tunable spring constant k 
connecting next nearest neighbor (NNN) sites [l6| and 
(2) increasing the volume fraction <p of pac ked spheres 
above the critical value </> c at jamming 
The isostatic lattices with their soft modes are then ap- 
proached continuously as k or A<p = (<f) — <f) c ) approach 
zero, and divergent length scales I*, vanishing frequen- 
cies uj* , and possibly vanishing shear moduli G (isotropic 
for jamming and the anisotropic modulus C44 = C xyxy 
for the square lattice as detailed below) can be identi- 
fied. In approach (2), the number of contacts increases 



asAz = z- z c ~ (A0) 1 / 2 , I* ~ (Az)-\ uj* ~ Az, and 
G ~ Az, whereas in approach (1) for the square lattice 



-1/2 



and G ~ k. 

In this paper, we investigate a third approach to iso- 
staticity in the square lattice: we populate NNN bonds 
with springs of spring constant K with probability V as 
shown in Fig. [TJ At nonzero V , the addition of an exten- 
sive number of NNN bonds removes all zero modes with 
a probability that approaches unity [20| as the number 
of sites A^ — s-oo, and as a result, the infinite lattice has a 
nonzero shear modulus for all V > 0. Thus, our model 
describes a rigidity percolation problem in which the per- 
colation threshold is at V = 0. It is the particular case 



2lL [22J of the more general rigidity percolation problem 
on a square lattice [23[ with AW and NNN bonds pop- 
ulated independently with respective probabilities Vnn 
and V in which T j nn = ^- This model shares underlying 
periodicity with approach (1) but it includes randomness 
analogous to approach (2). Adding a NNN spring in- 
creases the number of contacts by 1 so that V — [z—z^j/A, 
where z c = 4 in the AW square lattice. Unless otherwise 
stated in what follows, we use reduced units with k = 1 
and lattice constant a — 1 and unitlcss spring constants, 
clastic moduli, and frequencies: n/k— > k, Ga 2 /k — > G, 
and w/Vk — > uj. 

We study this random NNN model using the Coher- 
ent Potential Approximation (CPA) 23-25j], which gives 
good results for the conductivity of random networks 
near percolation p?j| and for rigidity percolation prob- 
lems [23[ except right in the vicinity of V = ? c , and we 
verify that it gives results that are in quantitative agree- 
ment with numerical simulations in our system. In the 
CPA, an effective medium of a uniform lattice with every 
NNN bond occupied by a spring with complex effective 
spring constant R m (u>) = K' m (ui) — ik'^uj), determined 
by a proper self-consistency condition, is used to capture 
the disorder average of the random lattice. From K m (ui), 
which is also equal to the complex shear modulus G(ui), 
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FIG. 1: (a) Square lattice with NN bonds with springs of 
spring constant k and NNN bonds with randomly placed 
springs with spring constant ft. The distortion depicted with 
dotted lines represents one of the zero modes of the lattice 
with no NNN springs, (b) Effective-medium lattice with 
springs of spring constant n m on all NNN bonds. In the 
CPA, the spring constant k 3 of a single NNN bond is changed 
to k or to with respective probabilities V and 1 — V ■ 

TABLE I: Dependence of I*, u>* , and G on V and Az. 

I* ~ p- 1 ~ (Az)- 1 lu*~V~Az G~P 2 ~(Az) 2 



we can calculate (following the procedures of approach 
(1) [lij]) the characteristic length I* and frequency uj* 
and the zero-frequency shear modulus G = R'{uj = 0), 
as summarized in Table |TJ As in the case of jamming, 
I* ~ 1/ui* ~ (Az) -1 , in agreement with the general cut- 
ting arguments of Ref. US]- The length /*, being the 
average distance between NNN bonds in any row or col- 
umn in the random lattice, marks the crossover from Id 
to 2d behavior in the effective medium, because NNN 
bonds couple neighboring Id rows or columns. The shear 
modulus, however, scales as G^V 2 ~ (Az) 2 , rather than 
as G ~ (Az) at jamming, implying highly nonaffine re- 
sponse near V = 0. If the response were affine, every 
equivalent NNN bond would distort the same way in re- 
sponse to shear, and G would be equal to Vk. Response 
becomes more nearly affine with G w Vk when ir 2 V ^> n. 
Figure [2] shows K m = G as a function of V for different k 
calculated from the CPA and via numerical simulations 
using the conjugate gradient method [28| to calculate the 
relaxed response of the system to an applied shear. 

The frequency dependence of k m (ui) is plotted in 
Fig. [3] In the nonaffine regime, it obeys a scaling law, 
k(lu) = K m h(uj/uj*), where h(w) approaches unity as 
w — > 0. k" {uj) vanishes as uj 2 at small w but be- 
comes nearly linear in uj for ui > 0.5a;*. This behavior 
corresponds to a shear viscosity that vanishes as uj at 
small u> but becomes a constant at large u>. A trans- 
verse phonon of frequency uj propagating along the y- 
direction (i.e., with q x = 0) has a wave number q{ui) = 
^/v^mM an d a mean-free path 1{uj) = y/ n' m (uj)T(uj), 
where r(w) = 2[K,' 1 ' n (ui)q 2 (uj)/uj}~ 1 is the decay time, im- 
plying that the Ioffe-Regel limit [13] q(uj)l{uj) = 1 occurs 
at 2k' m (uj) = k^(w), i.e., at uj w uj*. Thus uj* sets the fre- 
quency scale for the nearly isostatic modes and the scale 
at which plane-wave states become ill defined in agree- 




FIG. 2: (color online) Comparison of the CPA solution (lines) 
and numerical simulations on a 100 x 100 lattice (data points) 
for the effective medium spring constant K m as a function of V 
for k — 10 -2 , 10°, and 10 2 (in reduced units). Also shown are 
the nonaffine (n m = (tyV/2) 2 ) and affine limits (/t m = Vk). 
For the CPA at large V, we used the full dynamical matrix 
[Eq. {TJ] rather than the approximate forms of Eq. 




FIG. 3: (color online) Real and imaginary parts of h(uj/cj*) = 
b! — ih" (labeled respectively h' and h") and of K m (u)/c m for 
V = 10 -2 and 10 _1 (labeled respectively T, l", 2' and 2")for 
K = 1. Curves for V = 10" 3 and 10" 4 differ by less than 1% 
from the h curve and are not shown. The the full dynamical 
matrix [Eq. {l}] was used in the V = 10 _1 calculation. 

ment with recent studies of thermal conductivity near 
jamming [l9| . Because q y (uj*) ^ n/a, plane wave states 
with q x = are well-defined up to the zone edge. 

Because the zero modes on isostatic square lat- 
tice are uniform displacements of rows or columns, 
its phonon spectrum is identical to that of decou- 
pled one-dimensional chains with frequencies uj x ^ y (q) = 
2\smq x ^ y /2\ and density of states p(ui) = (2/n) /\fA — uj 2 
with a nonzero value l/ir at uj = as shown in Fig. SI 
When the effective-medium NNN coupling k m (uj) is 
added, the dynamical matrix becomes 

D xx {q) = D yy (q y7 q x ) = 4sin 2 (g 2 ./2)+4K m (w) sin 2 (g y /2) 
+ 4k to (w) sm 2 (q x /2) - 8R m (uj) sm 2 (q x /2) sm 2 (q y /2) 7 
D xy (q) = D yx (q) = 2k m (ui) s\i\{q x ) sm{q y ). (1) 

In the q — > limit, the dynamical matrix reduces to that 
of continuum elastic theory with D xx = Cuq 2 + Gnq 2 , 
where C\\ is a compression modulus and C44 the shear 
modulus. 644(01) is the complex shear storage modulus 
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FIG. 4: (color online) (a) Density of states p(oj) for (1) (green) 
a uniform lattice with k = K m on all NNN bonds, (2) (red) 
in the scaling nonaffine limit where k m (ui) — K, m h{uj /ui*), and 
(3) (Blue dotted) for V = 1CT 1 (4) (black dashed) Isostatic 2- 
mode limit of 2/tv w 0.64. (b) Density of states for a 100 x 100 
lattice with V = 10 _1 obtained via direct numerical calcula- 
tion (dots) and via the CPA (line) using the full rather than 
the approximate dynamical matrix of Eq. ((2J). Binning of the 
CPA result would wash out the spikes at low frequency at 
oj — q x = 27T7i/100 (for integer n). 



G(uj). Comparison of the continuum form with the small 
q limit of Eq. JT]) yields k m {uS) = G(uj). 

When |K m (oj)| <C 1, the off-diagonal terms in Dij can 
be ignored, and the low-frequency modes follow from 

D xx (q) f« q\ + ±k m {uj) sm 2 (q y /2) ps q\ + k m {uj)q 2 y (2) 

and a similar approximation for D yy (q). Replacing 
k m (u>) by its u> — > limit n m yields a characteristic 
length I* = \J\j 4k to through the comparison of q 2 with 
D xx (0, 7r) = An m and a characteristic frequency at the 

For q x > 1/1* 



zone edge of uj* = \J D xx (0, tt) = 2^J~k~n 
(or uj >uj*), the excitation spectrum is one-dimensional 
in q x . These observations along with K m ~ V 2 , which we 
derive below, lead to the results of Table HI 

To proceed with the CPA, we use the 2x2 phonon 
matrix Green's function of this effective medium 



G(q,o;) = [a; 2 I-D(q)] 



(3) 



In the CPA approximation [2J, [26j , an arbitrary NNN 
bond, say, between particles 1 and 2 as shown in 
Fig. [Tf b) , is replaced by a new one with a random spring 
constant n s with values n and with respective proba- 
bilities V and 1 — V. The dynamical matrix then changes 
to D v = D + V, where V is the potential given by [23J 



in real space, b = (e x + e y )/"j2 is the unit vector along 
the chosen NNN bond, and I and /' specify sites on the 
lattice. The potential V leads to a modification of the 
phonon Green's function, G^, (w), which can be calcu- 
lated following standard procedures: 

Gft, M = G ( _ ; , H + J2 G,_ {1 M ■ T luh ■ G la - V (u), (5) 

where Gi-i> is the Fourier transform with respect to q of 
G(q, u>) and where T = [1 — V-G] -1 • V is the scattering 
T-matrix. The effective spring constant k m (oj) is deter- 
mined within the CPA through the requirement that the 
average T vanish: V T| Ks=K + (1-V) T| Ks=0 = so that 

f{k m ,uj)k 2 m {io) - [1 + nf(k m ,ui)}k m (oj) + nP = 0. (6) 

The function / can be expressed as f(R m ,oj) = 
[2/(7t v / k^)]5(k„ 1 ,cj/ v /k^), where 



dq- 



1 _ e -Vrp(q,s) 



cos q 



(7) 



with p(q, s) = y^4sin 2 (g/2) — s 2 . In the limit r, s — >• 0, 
g(r, s) = 1, and thus f(k m , 0) -> [2/ (tt^/k^)} as n m -> 0. 
When y/rp(Tr,s) <C 1, the exponential in the numerator 
of <?(r, s) can be replaced by unity, and <?(0, s) = f;(s), 
g(s) 1 + (s 2 /8){ln[8/(Ves)] + «0/2)}. We expect K m 
to tend to zero with V so that in the small V limit, we 
can generally ignore the first term in Eq. ([6|). 

We consider first the static limit, ui = 0, for which the 
self-consistency equation for small V becomes 



-Tk = Q. 



The solution of this equation has two limits: 

if tt 2 "P < K, 
if 7T 2 -P > K, 




(8) 



(9) 



V 



,)(«i,i-*l,a)6®(^,i-*l',a)6, (4) 



as shown in Fig.[2J together with solutions of the full CPA 
equation (J6j) and numerical simulations. In the first 
case, K^n m 3> n m , and the solution for n m is obtained 
by ignoring the first term in Eq. ([5]): in the second case, 
the opposite is true, and K m is obtained by ignoring the 
second term in this equation. In the second case, every 
NNN bond distorts in the same way under stress, and 
response is affine. In the first (ttP/2) 2 < Vk, 

and response is nonaffine with local rearrangements in 
response to stress that lower the shear modulus to be- 
low its affine limit. Within the CPA, this result emerges 
because of the divergent elastic response encoded in G 
(and f(K m , 0)) as K m — > 0. As k approaches zero at fixed 
V , distortions produced by the extra bond decrease and 
the nonaffine regime becomes vanishingly small. 

For finite frequency uj. the effective medium spring 
constant is complex, k(ui) = k'(uj) — ik"(uj), where the 
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imaginary part k"(uj), which is odd in u and positive 
for bj > 0, describes damping of phonons in this ran- 
dom network. As in the static case, the nonaffine limit 
of the CPA result for at small V is the solution to 
k m f(k m ,uj) = V obtained from Eq. ((5]) by ignoring all 
but its last two terms. Following Eq. (|7|), at small k m and 
w, f(k mi uj) = [2/ (Tr^/k~^)]g(2 v / K m /k m u)/oj*). Thus in 
this limit, k m (uj) satisfies a scaling equation k m {uj) = 
K m h(uj/uj*). As oj -> 0, h(w) ->• 1 - w 2 {ln[4/(^ew)] + 
i(n/2)}, and k"(uj) ~ w 2 at small w. We calculated 
K m (u)/K m for P = 1CT 4 , HT 3 , 1CT 2 and HT 1 with the 
full CPA equation (j6)) and the nonaffine scaling function 
h(u>/cj*) for k = 1. The crossover from nonaffine to afhne 
behavior in the static limit is at "P = f/7r 2 w 10 _1 , so 
all cases but V = 10 _1 arc at or near the nonaffine limit. 
kJ[,(w) becomes greater than n' m {uS), and thus according 
the Ioffe-Regel criterion 27J, plane- wave phonon modes 
become heavily damped and ill-defined at ui « to* for all 
four values of V . 

The phonon density of states (DOS) p(uj), calculated 
from ImTrG Tn (q, uj) in the usual way, is plotted in 
Fig. @|a) as a function of lo/lu* . Curves for the three 
lowest V in Fig. SJa) collapse on to a common curve for 
u < 3u* . The curve for V = lO -1 departs from the com- 
mon curve at uj rj 0.5w* and is plotted in the figure. The 
large value of K," n {uJ*) in the random system removes the 
strong van Hove singularity at ui* of the uniform system. 
Figure SJb) compares the DOS for a finite lattice calcu- 
lated from CPA and by direct numerical diagonalization 
of the Hessian matrix using ARPACK [29| . The peaks in 
Figure BJb) at u) = q x = (2im/L) are due to finite size 
effects of the lattice with size L. 

Wc have used the CPA to analyze the static and dy- 
namic properties of a simple system on the threshold of 
isostaticity, namely a square lattice with NN springs and 
randomly distributed NNN springs. This system pro- 
vides clean analytic results about a random system near 
isostaticity, including nonaffine response near V = 0, 
and the scaling form for k m {ui) (which to our knowl- 
edge has not been observed in jamming systems), that 
can serve as a comparison point for more complicated 
systems. Our results strongly suggest that the divergent 
length I* ~ 1/uj* ~ (Az)" 1 is a common feature of all 
nearly isostatic systems in agreement with the arguments 
of Ref. They also unambiguously demonstrate that 
clastic moduli are not universal but depend on the geom- 
etry of the isostatic lattice. Further study is needed to 
determine exactly what properties of the isostatic lattice 
lead for example to a finite bulk modulus and a shear 
modulus vanishing as Az (as in jamming) or (Az) 2 (cur- 
rent system) or as (Az)° (kagome lattice j30|) or to one 



in which both B and G vanish as Az as in Ref. [31| . 
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